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Abstract 

Characterisations of the long-term behaviour of heart rate variability in humans 
have emerged in the last few years as promising candidates to became clinically sig- 
nificant tools. We present two different statistical analyses of long time recordings 
of the heart rate variation in the Eastern Oyster. The circulatory system of this 
marine mollusk has important anatomical and physiological dissimilitudes in com- 
parison to that of humans and it is exposed to dramatically different environmental 
influences. Our results resemble those previously obtained in humans. This suggests 
that in spite of the discrepancies, the mechanisms of long-term cardiac control on 
both systems share a common underlying dynamic. 
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1 Introduction 



The change with time in the size of the interval between two consecutive heart- 
beats is called heart rate variability (HRV). There are many (not necessarily 
independent) sources of HRV in people JjJ) (U), but the variations are largely 
controlled by the autonomic nervous system through the action of both the 
sympathetic and the parasympathetic branches, while the main mechanical 
influences are respiration and blood pressure. 



It has been found that long-term HRV shows 1/ f noise (|5|). This behaviour 
is found in many dynamical systems and draw attention to long-term HRV. 
The detrended fluctuation analysis (DFA) was introduced by Peng et al. Jfj) 
to study the long-range correlations found in HRV (7). A cross-over was iden- 
tified around a scale of 10 beats in all healthy subjects, signalling a change of 
dynamics when going from short to long time scales. Since then DFA has been 
used to characterise HRV in healthy conditions and in the presence of heart 
disease jsi jpl). Another promising approach is the analysis of HRV using 
wavelets, which is a mathematical technique specifically suited to analyse non- 
stationary series. The wavelet transform extracts the cumulative amplitudes 
of fluctuations of data at each point in time for a given scale (Q). Ivanov et 
al. (l3) presented the cumulative variation amplitude analysis (CVAA), where 
the inter-beat series were treated with consecutive wavelet and Hilbert trans- 
forms and an instantaneous amplitude is assigned to each inter-beat interval. 
It was found that the same Gamma distribution describes the distributions of 
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instantaneous amplitudes at all scales and for all the healthy subjects in the 
study. Further studies using wavelets on long-term recordings have explored 
the possibilit y t o define methods which could be used as markers of heart 



disease 



isiDiiit y t o 



The DFA and CVAA results suggest that there are intrinsic unknown dynam- 
ics underlying the long-term behaviour of the healthy human heart. It has 
been shown by Hausdorff and Peng (|2ll that it is extremely unlikely that 
the emergence of these complex patterns is due to having a system of many 
different independent influences each with their own timescale. The question 
of the origin of the universal long-term behaviour of HRV remains open. 

In this article we present DFA and CVAA studies of long-term HRV in the 
Eastern Oyster in conditions resembling those of their natural habitat ([22L In 
the first section, the circulatory system of the oyster is briefly described. The 
basic components of the system designed for the monitoring and acquisition of 
the heartbeat data are also presented. In the second section, the mathematical 
principles of DFA and CVAA are reviewed and then applied to the analysis 
of the oyster's heartbeat data. In the last section, the conclusions and some 
general remarks are given. 



2 The system under study 



2.1 The Eastern Oyster 



The Eastern Oyster is a fairly well studied mollusk (23$ |2J) which lives in 
coastal waters and lagoons from Canada to Mexico. It has an open circula- 
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tory system, i.e., the blood moves not only in the arteries and veins but also 
throughout the tissues. There are two accessory hearts which beat a couple 
of times per minute, independently from the principal heart. The main heart 
has three chambers. The two auricles receive the blood from the gills and send 
it, about half a second later, to the ventricle. The automatism of the heart 
is of a diffuse nature and contractions originate at any point of the ventricle. 
Contractions are not induced by impulses from the central nervous system 
and it is not known if there are any localised pacemakers. Two types of vesi- 
cles are found in the nerve endings of the myocardium, but it is not clear if 
they correspond to a rudimentary version of the sympathetic and parasym- 
pathetic systems ([25L The main external influences on the heart rhythm are 
the temperature, the level of oxygen and the salinity of water, while the main 
mechanical influences are the movement of the shell valves and the gills. When 
compared with people there is a factor of hundred in the size of the hearts, the 
period of respiration in the Eastern Oyster is two to three times longer and 
the heart beats about two times slower. The heartbeat receives perturbations 
from the shell valves and from the accessory hearts, which are not present in 
the case of the human heart. In summary all the suspected leading causes of 
the variability of the heart rate in people are not present or are quite different 
in the case of the Eastern Oyster. 



2.2 The measurements 



The experimental set-up was composed of a low-power laser diode, a fibre 
optic bundle, a photo-diode and the data acquisition and monitoring sys- 

nm. The radius of 
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the laser beam was 1 mm. The fibre optic bundle had a length of 1.3 m and 
a cross section of 38 mm 2 with a transmittance of approximately 60 percent 
at this wavelength. The active size of the EG&G Judson J16-5SP-R03M-SC 
germanium photo-diode was 3x3 mm 2 . The voltage signal was recorded with 
an AT-MIO 16 data acquisition card connected through a BNC-2080 multiple 
channel interface board. The system was controlled with a Lab View program 
written by us. A sketch of the technique is shown in Fig. 1. 

We have studied a set of six oysters with a uniform size of 7 cm and with a 
relatively thin section of the shell on top of the heart. Some of the oysters were 
measured several times, having as a result a set of 15 time series. There were at 
least 12 hours between measurements on the same mollusk. The oysters were 
kept in a big water tank under conditions resembling those of their natural 
habitat. A short time before the measurements, they were transfered to a small 
recipient on the focus of the optical set up. A system of pumps kept the water 
circulating and passing through other containers where the water was filtered 
and where the salinity, and temperature were controlled. 

3 The analyses 

The analyses shown below were performed on each of the inter-beat time 
series of the oyster's heartbeat. The laser light was pointed onto the beating 
heart. The intensity of the reflected light increases as the thickness of the walls 
increases in each systole of the ventricle. Hence, the inter-beat of a cardiac 
signal measured with our laser technique corresponds to the time between 
two consecutive peaks (see Fig. lb). The monitored signals of the oyster's 
heartbeat are: i) Highly non-stationary, ii) non periodic, and iii) irregular. 
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Further physiological perturbations on the cardiac signals such as gills and 
valves movements are clearly identified in the long-term. This set of special 
characteristics of our cardiac signals exclude the application of algorithms 
currently used to calculate the inter-beats in ECGs of humans. The most 



popular, are based on the QRS complex identification (j27l ). which does not 
apply to the cardiac signals from the oyster. 

To find the inter-beat intervals produced by the systole we have used the 
following method (see Fig. 2): (i) A set of boxes of the same width D is 
used to find the maximum in each interval, (ii) A secondary box of width d 
centred at the limit of two contiguous boxes is used as a range of confidence 
because the inter-beat period is not constant, (iii) The value of D and d are 
fixed visually such that the efficiency for finding the peaks is optimal, (iv) 
Furthermore, after the calculation of the inter-beat intervals, the outliers are 
filtered with a similar procedure as suggested in Ref. 



28l ). We obtained as 
a result of applying the algorithm, 15 time series with lengths varying from 
10403 to 25802 inter-beats. 



3.1 Detrended Fluctuation Analysis 

The detrended fluctuation analysis is a technique that permits to identify long- 
range correlations in non-stationary time series. As commented in Ref. (29), 
DFA has been applied to the study of a broad range of systems, such as the 
human gait, DNA sequences, the heartbeat dynamics, the weather, and even 
in economics. Specially, in the analysis of natural inter-beat heartbeat fluctu- 
ation, DFA has helped to discriminate healthy from heart diseased humans. 
It has provided also, a quantitative difference between old and young people. 
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Several works show the robustness of DFA, although improvements are still 
being done DFA is a simple yet powerful tool for studying physiological 
data. 



a. 



The DFA is as follows (p). Let x(i) be a time series. Then (i) integrate x(i), 
(ii) divide the time series in n equal amplitude intervals, (iii) in each box 
of width n do a polynomial fit of order / (it defines DFA-/ analysis): s l n (k), 
(iv) eliminate the polynomial trend in each box, (v) calculate the root mean 
squared fluctuation F as a function of the n intervals, and (vi) do steps (i)-(v) 
for several box widths to find the functional relation between F(n) and n. 

The presence of scaling in the original signal produces a straight line in a 
double log plot of F{n) versus n. If the slope, a, is 0.5, the data is uncorrected 
and corresponds to a random walk jsii. A slope between 0.5 and 1 signals 
the presence of a long-range power law, where a = 1 corresponds to 1/f 
noise ([32L For a slope bigger than 1, the correlations no longer correspond 
to a power law. A value of 1.5 indicates Brownian noise (|33). Healthy people 
have a slope of 1.5 for values of n lessjthan 10 beats, and a slope of 1 for time 
scales between 100 and 10000 beats 

Using DFA-1 we found (Fig. 3) that all the oysters in our study present 1/f 
noise behaviour for scales above log(n) ~ 2.5, corresponding to ~ 300 beats. 
The average slope in this region is «2 = 1-08 ± 0.14 where the error is the 
standard deviation of all the samples. For shorter time scales, between 10 and 
100 beats, a slope a± = 0.61 ± 0.04 is found. The slope at even shorter scales 



(log(n) < 1) is also close to 1 in all cases, but the D 



tially large intrinsic systematic effects in this range (|29l). making it difficult to 



A method has poten- 



extract reliable information. We found the same variation in the results from 
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measurement to measurement when comparing different sets from the same 
oyster and sets of data from different oysters (see Table 1). It is interesting to 
note that the respiration and the accessory hearts have independent oscilla- 
tions with periods around 20 to 40 beats. In contrast, the shell valves present 
intermittent activity but, when active, the period is also in the 20-30 beats 
range. These complex perturbations on the oyster's heartbeat could be the 
source of the aci slope at short time scales (34; 3^). To discard the possibility 
that polynomial trends could be the source of the cross-over (3^) obtained 
in this work, we performed analyses of oyster's heartbeat using DFA-Z for 
/ = 2,3,4, and found in each case a cross-over approximately in the same 
region indicated by DFA-1, although cti decreased down to 20% while ct\ 
keeping its same value. 



3.2 Cumulative Variation Amplitude Analysis 

We also analysed our measurements using CVAA. This technique that is based 
on consecutive wavelet and Hilbert transforms was applied for the first time 
in the study of natural heartbeat fluctuation. It was found that a common 
Gamma distribution characterises a group of healthy people. On the other 
hand, in the case of a group of people suffering sleep apnea, it was found that 
such data collapse does not happens. Even more, in some cases it was not 
possible to get a Gamma distribution (ljj). In general, a Gamma distribution 
is characteristic of physical systems out of equilibrium. Hence, the results 
previously commented suggest that the heartbeat in healthy people owns an 
intrinsic underlying dynamics. 

Mathematically, CVAA consists of the next steps: (i) Choose adequate scales 
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to analyse the data, (ii) from the original series, a set of /eries each at a 
different scale is obtained using a continuous wavelet transform. There are 
many wavelet families to choose to perform this step and several have been 
tried. Each family eliminates local polynomial trends from the signal in a 
different way. The coefficients c of the wavelet transformation in each scale 
reflect the cumulative variation of the signal, (iii) Then, each of the new time 
series is processed with a Hilbert transform to extract the instantaneous am- 
plitudes h of the variations at each point in the series, (iv) Construct the 
time series y = c + ih and calculate the amplitudes A = ^/c 2 + h 2 , (v) fi- 
nally, the histogram of these amplitudes is normalised to 1 to form a prob- 
ability distribution, P(x), which is then re-scaled such that x —>■ xP max and 
P{x) -> P{x)/P max . 

Remarkably, we found that each distribution of instantaneous amplitudes is 
fitted by a Gamma distribution jisi (Fig. 4a). Furthermore, as in the case 
of healthy people, the distributions for all the oysters in the study are well 
described by the same Gamma distribution (Fig. 4b), i.e., there is a common 
parameter v J^) which describes the normalised distribution of instantaneous 
amplitudes from any oyster. This behaviour is found at scales 2 3 ' 4 ' 5,6 and for all 
the wavelets analysed Daubechies (moments 3-10), Gaussian (moments 
3-10), Meyer, Morlet, and B-spline Biorthogonals (decomposition moments 1 
and 3). The results obtained for the parameters of the Gamma distributions 
were all very similar. As an example of the robustness of our results, in Table 2 
are shown the values of the parameters of the fits to the Gamma distributions 
performed with orthogonal, biorthogonal, and non-orthogonal wavelets. 

The numerical value for the v parameter in the Eastern Oyster is v = 1.0 ±0.2 
which is slightly lower than the value of v = 1.4 ±0.1 found for healthy people 
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during sleep hours by Ivanov et al. using Gaussian wavelets (15 



4 Conclusions 

Using the DFA and CVAA methods we find long-range correlations and scaling 
in the long-term HRV of the Eastern Oyster. DFA shows 1/ / noise behaviour 
at large scales and a cross-over to a smaller slope for scales of the order of 300 
beats for all oysters in the study. The cross-over happens at a scale well above 
the region where DFA presents some bias. The source of the cross-over seems 
to be the result of the complex interactions between the components of the 
circulatory system of the Eastern Oyster. Models of this phenomena such as 
polynomial trends added linearly to a correlated signal do not seem to apply 
to the case of the Eastern Oyster. 

With CVAA we find that all oyster records collapse to a Gamma distribution 
with the same numerical value of the v parameter, v w 1, for a wide variety 
of wavelets and scales. These results are remarkably similar to those previ- 
ously reported in the study of healthy people, in spite of the fact that the 
circulatory system of the Eastern Oyster and the influences it is exposed to 
are dramatically different from those in the case of people, pointing thus to 
an intrinsic origin of these complex patterns. Characterisations of long-term 
HRV are promising candidates for clinical prognostic tools making it vital to 
understand its origin in order to exploit fully this type of techniques. Our 
results pose stringent constrains and offer new hints and challenges to models 
attempting to describe the long-term dynamics of the heart. 



10 



Acknowledgements 

We thank D. Vera and J. Bante for their technical assistance and G. Oskam 
for fruitful discussions. This work was partially supported by Conacyt Grant 
28387E. 

References 

[I] Task Force of the European Society of Cardiology and the North American 
Society of Pacing Electrophyisiology, M. Malik et al. Heart rate variability: 
Standards of measurement, physiological interpretation and, clinical use. 
Circulation 93 (1996) 1043. 

[2] D. T. Kaplan, M. Talajic, Chaos 1 (1991) 251. 
[3] S. Akselrod et al, Science 213, (1981) 220. 

[4] J. J. Goldberger, S. Challapalli, R. Tung, M. A. Parker, A. H. Kadish, 

Circulation 103 (2001) 1977. 
[5] M. Kobayashi, T. Musha, IEEE Trans. Biomed. Eng. BME 29 (1982) 456. 
[6] C. -K. Peng, S. Havlin, H. E. Stanley, A. L. Goldberger, Chaos 5 (1995) 

82. 

[7] C. -K. Peng et al., Phys. Rev. Lett. 70 (1993) 1343. 
[8] P. Ch. Ivanov et al., Europhys. Lett. 48 (1999) 594. 
[9] A. Bunde et al., Phys. Rev. Lett. 85 (2000) 3736. 

[10] S. M. Pikkujamsa, T. H. Makikallio, K. E. J. Airaksinen, H. V. Huikuri, 
Am. J. Physiol. Heart Circ. Physiol. 280 (2001) H1400. 

[II] Y. Ashkenazy et al., Phys. Rev. Lett. 86 (2001) 1900. 
[12] J. W. Kantelhardt et al., Phys Rev E 65 (2002) 051908. 
[13] J. C. Echeverria et al., Chaos 13 (2003) 467. 



11 



[14] I. Daubechies, Ten Lectures on Wavelets, Society for Industrial and Ap- 
plied Mathematics, Philadelphia, PA, 1992. 
[15] P. Ch. Ivanov et al, Nature 383 (1996) 323. 

[16] S. Thurner, M. C. Feuerstein, M. C. Teich, Phys. Rev. Lett. 80 (1998) 
1544. 

[17] L. A. Nunes Amaral, A. L. Goldberger, P. Ch. Ivanov, H. E. Stanley, 

Phys. Rev. Lett. 81 (1998) 2388. 
[18] S. Thurner, M. C. Feuerstein, S. B. Lowen, M. C. Teich, Phys. Rev. Lett. 

81 (1998) 5688. 
[19] V. Pichot et al, J. Appl. Physiol. 86 (1999) 1081. 

[20] G. McCaffery, T. M. Griffith, K. Naka, M. P. Frennaux, C. C. Matthai, 

Phys. Rev. E 65 (2002) 022901. 
[21] J. M. Hausdorff, C. -K. Peng, Phys. Rev. E 54 (1996) 2154. 
[22] P. A. Ritto, Ph. D. thesis, Centro de investigation y de estudios avanzados 

del IPN, Merida, Mexico, 2003. Unpublished. 
[23] P. S. Galtsoff, The American oyster Crassostrea virginica Gmelin, Fishery 

Bulletin of the Fish and Wildlife Service, Washington, 1964. 
[24] V. S. Kennedy, R. I. E. Newell, F. Eble Albert Eds., The Eastern Oyster 

Crassostrea virginica, Maryland Sea Grant College, College Park, 1996. 
[25] P. G. Beninger, M. Le Pennec, in Scallops: Biology, Ecology Aquaculture. 

S. E. Schumway Ed., Developments in Aquaculture and Fisheries Science, 

Volume 21, Elsevier, 1991 pp. 133. 
[26] P. A. Ritto, J. G. Contreras, J. J. Alvarado-Gil, Meas. Sci. Technol. 14 

(2003) 317. 

[27] A. L. Goldberger et al., Circulation 101 (2000) 215, and references 
therein. 

[28] K. Ho et al., Circulation 96 (1997) 842. 



12 



[29] J. W. Kantelhardt, E. Koscielny-Bunde, H. H. A. Rego, S. Havlin, A. 

Bunde, Physica A 295 (2001) 441. 
[30] H. Yang et al., preprint |c7md-mat/ 0201206] (2002). 

[31] E. W. Montroll and M. F. Shlesinger, in Nonequilibrium Phenomena II. 

From Stochastics to Hydrodynamics, edited by J. L. Lebowitz and E. W. 

Montroll, North-Holland Amsterdam, 1984, pp. 1-121. 
[32] P. Bank, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381. 
[33] C. -K. Peng, S. Buldyrev, A. L. Goldberger, S. Havlin, F. Sciortino, M. 

Simons, and H. E. Stanley, Nature 356, (1992) 168. 
[34] Z. Chen, P. Ch. Ivanov, K. Hu, H. E. Stanley, Phys. Rev. E 65 (2002) 

041107. 

[35] K. Hu et al, Phys. Rev E 64 (2001) 011114. 

[36] Let x be a real variable. The Gamma distribution function is defined as 
P u (x, u) = u u+1 x u exp(— ux)/T(v + 1) where v is a real positive integer and 
u = u/xo, where xq localises the peak of the distribution. 
[37] The Mathworks, Inc., Wavelet Toolbox, WWW page 
:/ /www.mathworks.com/access/helpdesk/help/toolbox/wavelet , 
2004. 



13 



FIGURE CAPTIONS 



Figure 1. (a) A low power laser is pointed to the beating heart. The walls 
reflect more light when contracted than during the diastole, (b) The periodic 
variation of light intensity is measured with a photo-diode whose voltage 
output varies with time capturing the beating of the oyster's heart J3). Here 
the stars show the peaks found by our algorithm. 

Figure 2. This is a schematic representation of the way our algorithm for 
finding heartbeat peaks works. The cardiac signal is divided with a set of 
boxes of equal length D. Due to the fact the heartbeat period is not constant, 
a secondary set of boxes of width d < D are localised at the limit of two 
contiguous boxes. In the plot D = 1.25 and d = 0.5 are common values which 
optimises the search of peaks. See text for more details. 

Figure 3. (a) Result of the DFA performed on an inter-beat series from an 
oyster showing the cross over behaviour. The error bars represents the fluc- 
tuations in F(n) in a region ±0.05 around each n. The slope for longer time 
scales, «2 corresponds to the 1/ / noise behaviour which is also found in healthy 
human hearts. For shorter scales a behaviour close to white noise is found in- 
dicating that the signal is almost completely random in these time scales, (b) 
The value of aii and ct2 for all the samples. The solid lines correspond to a± 
and «2. The error bars reflect the variation of the slope when changing the 
start and end points of the fitting range. The bigger error bars in reflect 
the statistical fluctuation at large n. Note that the cross over does not happen 
in a point but that there is a transition region (Table 1). 

Figure 4. (a) Result of the cumulative variation amplitude analysis performed 
on a heart series from an oyster at a scale corresponding to 2 4 — 1 beats using 
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the fourth wavelet of the Daubechies family. The points are the data, while 
the solid line is the result of a fit to a Gamma distribution. The parameters 
obtained from the fit are v = 1.02 ± 0.01 and x = 0.37 ± 0.01 with a x 2 /dof 
of 0.3. (b) The data points, corresponding to the same scale and wavelet, for 
all oyster records. All of them collapse to a single Gamma distribution. The 
solid line corresponds to the same value of the parameters v and x shown in 
a (Table 2). 
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FIGURE 4b 
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Table 1: Cross over region and slopes before (a±) and after (02) the cross over for 
each sample. The error bars reflect the variation of the slope when changing the 
start and end points of the fitting range. The bigger error bars in 0.2 reflect the 



fluctuation at large n. 


Sample 


n of cross over region 






1 


2.2-2.7 


0.58 ± 0.02 


1.05 ± 0.05 


2 


2.2-2.7 


0.58 ± 0.01 


1.10 ± 0.15 


3 


2.4-2.9 


0.58 ± 0.01 


1.20 ± 0.15 


4 


2.1-2.4 


0.53 ± 0.04 


0.95 ± 0.10 


5 


2.5-2.8 


0.58 ± 0.03 


0.90 ± 0.10 


6 


2.4-2.6 


0.60 ± 0.02 


0.95 ± 0.10 


7 


2.0-2.2 


0.64 ± 0.01 


1.10 ± 0.10 


8 


2.4-2.8 


0.59 ± 0.01 


0.95 ± 0.10 


9 


2.4-2.7 


0.66 ± 0.02 


1.30 ± 0.15 


10 


1.9-2.1 


0.67 ± 0.02 


1.15 ± 0.10 


11 


2.2-2.5 


0.63 ± 0.01 


1.05 ± 0.05 


12 


2.1-2.6 


0.75 ± 0.01 


1.15 ± 0.20 


13 


2.3-2.8 


0.66 ± 0.06 


1.25 ± 0.10 


14 


2.3-2.8 


0.63 ± 0.05 


1.25 ± 0.10 


15 


2.1-2.5 


0.63 ± 0.03 


0.85 ± 0.15 
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Table 2: Average values for the parameters of the Gamma distribution calculated 
for each wavelet at each scale from the fitted values to each separated sample. The 
quoted errors are the sum in quadrature of the individual errors as given by the fit. 
Note that in this sample are included results using three types of wavelets: Non- 
orthogonal (Gaussian-3), orthogonal (Meyer and Daubechies-4) , and biorthogonal 
(B-spline-1.3. Nomenclature used is as follows J37). Right (left) index corresponds 
to decomposition (reconstruction) moment). 



Wavelet 


Scale 


V 


Xq 


Wavelet 


Scale 


u 


Xq 


gaus3 


3 


0.94 ± 0.36 


0.36 ± 0.07 


biorl.3 


3 


1.11 ± 0.56 


0.42 ± 0.25 


gaus3 


4 


0.94 ± 0.35 


0.36 ± 0.08 


biorl.3 


4 


0.94 ± 0.54 


0.36 ± 0.10 


gaus3 


5 


0.91 ± 0.36 


0.35 ± 0.08 


biorl.3 


5 


0.98 ± 0.80 


0.37 ± 0.12 


gaus3 


6 


0.89 ± 0.36 


0.35 ± 0.09 


biorl.3 


6 


0.93 ± 0.40 


0.35 ± 0.08 


meyr 


3 


1.09 ± 0.32 


0.39 ± 0.05 


db4 


3 


1.01 ± 0.38 


0.38 ± 0.07 


meyr 


4 


1.07 ± 0.28 


0.39 ± 0.05 


db4 


4 


1.01 ± 0.29 


0.37 ± 0.06 


meyr 


5 


1.06 ± 0.33 


0.39 ± 0.06 


db4 


5 


0.99 ± 0.33 


0.37 ± 0.06 


meyr 


6 


1.01 ± 0.36 


0.37 ± 0.07 


db4 


6 


0.97 ± 0.34 


0.37 ± 0.07 
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